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Recent numerical calculation of the intrinsic thermal Hall conductivity of nodal d-wave supercon¬ 
ductors in the mixed state revealed a rapid increase of this quantity above an onset temperature. 
Interestingly, this defines a measurable energy scale in an otherwise gapless state. Using the mathe¬ 
matics of magnetic coherent states, in this paper such energy scale is related to a dynamical process 
associated with the Andreev scattering of an electron wavepacket moving along the constant energy 
contours in the momentum space. This energy scale is then used to obtain an improved scaling col¬ 
lapse of numerically calculated thermal Hall conductivity in a tight-binding model as a function of 
temperature, magnetic field and the d-wave pairing amplitude at various band fillings. The results 
indicate that the mentioned onset temperature is associated with the ability of the quasiparticle 
wavepacket to complete its semiclassical orbit before it is appreciably scattered by the supercon¬ 
ducting condensate. 


I. INTRODUCTION 


The electrical Hall effect is an important technique in 
materials characterization. Unfortunately, it provides lit¬ 
tle useful information below the superconducting transi¬ 
tion temperature, even in type II superconductors, for 
which the magnetic field penetrates the bulk of the sam¬ 
ple. This is because no transverse voltage can be estab¬ 
lished in a superconductor, assuming, as is done thought 
this paper, that it superconducts i.e. that the vortices are 
rigidly pinned and not driven into the flux flow regime^. 


On the other hand, a superconducting sample in an 
external magnetic field H and subject to a small heat 
current density jQ, may exhibit a thermal Hall effect, 
i.e. a temperature gradient perpendicular to both H and 
jQ. The thermal Hall conductivity, K xy , is then defined 
as 3Q X ~ ~ K ’ x yl[f- 1 11 the case of extreme type II su¬ 
perconductors considered here, the magnetic field inside 
the sample is practically uniform, but, because the ele¬ 
mentary (Bogoliubov) quasiparticle excitations inside a 
superconductor are a coherent superposition of an elec¬ 
tron and a hole^, the usual theory of thermal Hall effect 
in normal metals^ does not apply directly. Development 
of such theory is therefore an important step towards ex¬ 
tending Hall measurements into the realm of supercon¬ 
ductivity. 


In a model of non-interacting Bogoliubov quasiparti¬ 
cles, the intrinsic contribution to K xy can be related to 
the energy dependence of the quasiparticle current Hall 
response^. The intrinsic contribution is the part inde¬ 
pendent of the impurity scattering; it is finite and well 
defined without any impurities and is expected to domi¬ 
nate in the clean limit. In the superconductor, the quasi¬ 
particle current is distinct from the electrical current^: if 
the quasiparticle Hamiltonian operator is H , the former 


is proportional to the quasiparticle velocity 


r ,H 


/ih, 


and the latter to dH/d A. Moreover, if the vortices are 


arranged in a perfect lattice, then Bloch theorem can be 
employed^, and the thermal Hall conductivity can be re¬ 
lated to the energy dependence of the Berry curvature of 
the quasiparticle sub-bands in the (vortex) crystal mo¬ 
mentum Brilluoin zone^. Using such approach, it was 
recently shown that at low magnetic field H , the intrin¬ 
sic contribution to n xy exhibits a simple scaling with H, 
and shows a rapid increase from negligible values at low 
temperature to values of order l/H at a characteristic 
onset temperature^. In the model used in Refill, the on¬ 
set temperature was shown to increase with increasing 
A, the pairing amplitude of the tight-binding lattice d- 
wave superconductor whose H = 0 pairing function is 
2 A (cos — cos k y ). While this successfully captures the 
most important dependence of the onset temperature of 
K xy , the numerical results of the Ref[H at fixed band fill¬ 
ing displayed additional (weak) dependence on the Dirac 
cone anisotropy (see Fig. 2 of Ref]E). This feature has 
not been explained. In addition, as found in this work, 
there is an additional dependence of the onset tempera¬ 
ture on the band filling in the tight-binding model used 
in the Refill. 

As explained in this paper, such features are a conse¬ 
quence of the particular lattice model adopted in Refjj|; 
the onset temperature dependence on A reported therein 
indeed captures the main essence of the effect. The men¬ 
tioned residual dependence can be naturally understood 
by picturing high energy quasiparticle wavepackets semi- 
classically moving along the contours of constant (nor¬ 
mal) energy. As the energy of the quasiparticle is low¬ 
ered, the amplitude of its Andreev scattering increases 
along the anti-nodal portions of its contour, and at some 
point becomes prohibitively large for the wavepacket to 
complete its semiclassical orbit before it is apprecia¬ 
bly scattered by the superconducting condensate. This 
marks the energy scale e*, which obviously increases with 
increasing A. However, due to the tight-binding disper¬ 
sion and pairing function used in the model of RefJfj, e* 
has additional dependence on the Dirac cone anisotropy 




as well as the band filling. To illustrate the band fill¬ 
ing dependence, consider the magnitude of the pairing 
function on the Fermi surface in the anti-nodal direction, 
2A(1 — cos fcf); its value depends not only on A, but also 
on ItF which depends on the band filling. If, instead of 
using A to rescale the temperature, e* is used, then the 
family of curves for a range of values of the Dirac cone 
anisotropy and band Filings collapses onto a single scal¬ 
ing curve (see Fig. HJ). Such improved scaling - combined 
with the explicit calculation for the scattering amplitude 
formulated in continuum and using magnetic coherent 
states presented below - therefore strongly supports the 
above physical picture. It also indicates that n xy may 
be a way to measure the ability of the quasiparticles to 
complete their semiclasical orbits before they are appre¬ 
ciably Andreev scattered, providing useful spectroscopic 
information about unconventional superconductors. 

The primary focus of this paper, just as in RefjB, is 
the limit Hco c <C A <C Ep, where the Fermi energy Ep 
is to be measured from the band minimum or maximum, 
whichever gives the smaller value, and w c = eH /me is the 
cyclotron frequency of a point particle with charge e and 
mass m. In this regime a naive perturbation theory in 
A would appear to break down. However, as mentioned, 
the key insight advanced here is that the high energy 
states must be weakly affected by the pairing term, de¬ 
spite their separation - set by Hu> c - being much smaller 
than the pairing term amplitude. In order to obtain the 
energy scale where A ceases acting perturbatively, a first 
order time dependent perturbation theory calculation is 
performed using as the starting state a magnetic coherent 
stated. Such states are exact eigenstates of the time de¬ 
pendent Schrodinger equation in a uniform magnetic field 
in the symmetric gauge, but they are not the stationary 
states - starting with a stationary state is often assumed 
in the quantum mechanics textbooks explaining time de¬ 
pendent perturbation theory, but it is, of course, not 
necessary^. A magnetic coherent state describes a Gaus¬ 
sian wavepacket moving along circular trajectory in the 
real space with the angular frequency ui c and the width of 
the Gaussian set by the magnetic length ip = y/hc/eH. 
In the absence of any other perturbations, the wavepacket 
width does not change in time. If one writes the Hamilto¬ 
nian operator for the electron, (p— |A) 2 /2m, in terms of 
the harmonic oscillator ladder operators as hui c (a^a+ ^), 
and that of the hole, (p + |A) 2 /2to, as hu c {b^b + |), 
then the magnetic coherent states | a/3) are the simulta¬ 
neous eigenstates of a and b. At finite time the solution 
of the time dependent Schodinger equation is |e _ “ c *o;, j3) 
when the dynamics is generated by (p — -A) 2 /2m, and 
\a, e~ lu>ct /3) when by (p T |A) 2 /2m. 

Thus, the main finding presented in this paper is that, 
as long as fuv c A <C Ep - where on the tight-binding 
lattice with a unit lattice spacing and with the hopping 
amplitude t, Huj c should be understood as - and, 

as long as the filling does not coincide with the vicinity of 
the van Hove singularity, the thermal Hall conductivity 


has the scaling form 


Kxy(T, A, H, n) 


= k xv (T,0,H,h) x T 



where n xy (T,0, H, fi) is the clean limit normal state 
thermal Hall conductivity, which obeys free Fermion 
Wiedemann-Franz law, and scales as ~ T/H. Here, kp is 
the Boltzman constant, which will be set to unity in what 
follows, unless stated explicitly otherwise. The energy 
scale e* depends on the magnetic field H only through 
a possible 17-dependence of A and /i. e» is to be deter¬ 
mined as follows: consider the normal state dispersion 
and the pairing function Ak- Then, as we move along 
the closed contours of constant |ek — ji\ shown in FigjTl 


the quantity 


Ak 


measuring the amplitude of Andreev 


scattering, varies. For the d-wave superconductor of in¬ 
terest here, this quantity is peaked in the antinodal di¬ 
rection. As we approach the Fermi level, there are two 
contours of constant |ek — /r|, one inside and one outside 

the Fermi surface, for which the peak value of is 


equal to a pure number 0* of order unity which will be 
specified shortly. Then, as shown in Fig. |T| e* is the 
lesser of the two such values of |ek — /i|. 


For the specific case of lattice d-wave superconduc¬ 
tor considered here, at H = 0 the pairing amplitude 
is Ak = 2A(cosfe r — cos k y ) and the normal state dis¬ 
persion is £k = — 2t(cosfca; + cosky). This results in 

e» = A ■ As shown in Fig[2l the scaling collapse 

of K xy is achieved for 0* ss 0.4, a value which is inter¬ 
estingly close to l/\/27r. The resulting scaling function 
J-(x ) is monotonically increasing, approaches 1 for large 
x , and displays a rapid onset at x ~ 0.1 (see Fig. [2]). It 
should also be mentioned that, because we are interested 
in the limit Hlo c -C A -C Ep, the Zeeman effect, corre¬ 
sponding to a trivial shift of all quasiparticle energies, is 
ignored here. 


The rest of the paper provides details of the calcula¬ 
tions which lead to the above assertions. In Section II, 
the mathematics of the magnetic coherent states in sym¬ 
metric gauge is reviewed. The methods for construct¬ 
ing the pairing order parameter in the vortex state, and 
in the symmetric gauge, are reviewed in Sec Ha. The 
time dependent perturbation theory to first order in the 
pairing term, and in the basis of the magnetic coherent 
states, is described in Sec lib. The entire formulation 
in Sec II is in continuum. The formulation on the dis¬ 
crete tight-binding lattice, along with the formula used to 
numerically compute the thermal Hall conductivity from 
numerically diagonalizing the tight-binding Hamiltonian, 
are reviewed in Sec III. Discussion is in Sec IV, and the 
details of the perturbative calculation with the magnetic 
coherent states are delegated to the Appendix. 












FIG. 1. Illustration of the physical process which determines 
the energy scale, e», for the onset of the intrinsic thermal Hall 
conductivity K xy in the lattice d-wave superconductor. In this 
figure, representing the 1 st Brillouin zone, the pairing ampli¬ 
tude is taken to be Ak = 2 A (cos k x — cos k y ), and the normal 
state dispersion is ek = — 2t(cos k x +cos k y ). The shaded con¬ 
centric contours are the contours of constant |ek — /t|, with 
fj, = 2t. The Fermi surface (FS), where |ek — fJ,\ = 0 is 
marked (yellow). The two d-wave shaped lines (red) corre¬ 
spond to contours of constant |Ak/(ek — yu) | = #„. The value 
of 6* = -j== S3 0.4 chosen here is the same as in Fig. [2] where 
it is shown to result in the collapse of the numerical data. 


II. MAGNETIC COHERENT STATES 


This section follows the original article by Malkin and 
Man’ko^. It is included here in order to establish notation 
and the main mathematical identities which will be used 
in later sections. This formulation is in continuum. 

In the symmetric gauge A = \Hz x r = |H (—y, x, 0). 
The cyclotron frequency is uj c = eH / (me), and let t = 
\/hc/eH. Note that this differs by a factor of 1/ y/2ii from 
the definition of the magnetic length £h = y^hc/eH used 
earlier. 

The Schrodinger Hamiltonian operator for the electron 
in the magnetic field is then 


n = 



2 m 


( 2 ) 
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FIG. 2. Scaling of the thermal Hall conductivity in the 
mixed state of the lattice d-wave superconductor K xy = 
n xy (T, A, H, fi), re-scaled by its (free Fermion) normal state 
value K'xy = K xy (T, 0, H, fj,). The temperature is rescaled by 
e„, an energy scale associated with Andreev scattering, as 
discussed in the text. Various values of the Dirac cone dis¬ 
persion anisotropy a = vf/va = t/A and chemical potential 
/r are shown in the legend. The magnetic length £h = 28a in 
all these calculations, with the square vortex lattice oriented 
along tight binding unit cell diagonal (see inset of the Fig 2 
of Ref;5j where magnetic field scaling has been established). 


1, where 


Then, 


V2 


e + 


dC 


^ = -4 


V2 




H = hjJc 




( 4 ) 

( 5 ) 

( 6 ) 


Note that there is another set of raising and lowering 
operators, satisfying [b,^] = 1, where 


b 

6 t 


1 

i 

V2 


f + 


e- 


dt 

d 

d£* 


( 7 ) 

( 8 ) 


These do not appear explicitly in the Hamiltonian, but, 
importantly, they commute with the previous ones: 

[a, b] = [a, b f ] = 0. (9) 


Let us define a dimensionless variable 

€=(x + iy)/( 2£), (3) 

and raising a and lowering operators satisfying [a, a^] = 


They therefore represent a constant of motion. For a par¬ 
ticle with opposite charge, the Schrodinger Hamiltonian 
operator can be written as 




2 


2m 



= hjJ, 


( 10 ) 

















The ’’vacuum” state 100) is simultaneously annihilated 
by a and 6, and in coordinate representation is given by 

<r|00) = (11) 

This state is used to build coherent states^. In order to 

do so, define the unitary operators 

D(a) = e Qat - Q * a , (12) 

D(/3) = (13) 


Note that the exponents commute. This operator com¬ 
mutes with (p — —A) , whose ground state wavefunc- 
tion 

—^ 2 =— e -(* 2 +J/ 2 )/( 4 ^) (20) 

serves to generate the order parameter; more precisely 
and as discussed below, its center-of-mass coordinate de¬ 
pendence. 

The irreducible representation for the magnetic trans¬ 
lation group (see e.g. Ref. EH) are 


where a and /3 are two complex c-numbers. Clearly, the 
two operators commute: 


D(a),D((3) 


= 0. 


(14) 


The common coherent state of a and b is 

\ap)=D(a)D(f3)\00). (15) 


In the coordinate representation, such state has the form 

(rla/3) = -i^e -c * f e^ t+iV5ar e _<a/J e _ ^( |a|2+lfl|2 ). 
x 1 ' y/toil 

(16) 

This is a Gaussian centered at £ = (/3* + ia) / y/2 
and modulated by the phase which grows linearly with 
x and y. The kinetic momentum for the electron, 
p x — | A x + i (p y — | Ay) = y/2ha/i, in such a state is 
peaked at \/2ha/£. The kinetic momentum for the hole, 
p x + -A x + i (p y + |A y ) = y/2ihtf /£, in such a state is 
peaked at ih/3* /£. 

The coherent states form an overcomplete set, and can 
be used to construct the resolution of identity 

f da* da f df}*d/3 , , 

J 97r - J -(r\a(3){a(3\r ) = <5(r- r ). (17) 

Here f ( \ _ f°° f°° dtftea d^sma / \ 

J 2ni '•***■' J—oo J—oo 7r 


A. Pairing order parameter in symmetric gauge 

The development in this section follows the work of 
T. Kita— . Because w e are de aling with charge 2e order 
parameter, let £» = y/hc/2eH = £/y/2. Then, consider a 
set of 2D lattice points 

R = mai -(- n 2 a 2 , (18) 

where n\ and 712 are integers. The primitive lattice vec¬ 
tors are ai = (ai x ,ai y ,0) and a 2 = ( 0 , 02 , 0 ), where 
ai x a 2 = 27r£ 2 . 

In the symmetric gauge, the operator which corre¬ 
sponds to the translation by a lattice vector R, followed 
by a gauge transformation, is 

T(R) = e -i(fl»*--R*i/)/(2^) e - R -v_ ( 10 ) 


L>( q )(R) = m 2 (2!) 


Then, at q = 0, the s-wave Abrikosov order parameter 
can be written as 


A^e" nin2 T(R), 


rM* 2+ y 2 ) 


R 


oo oo 


A E £ < 

n 1= —oo n 2 — — oo 


( 22 ) 


where 

j. H - 'iRy 

CR = V2L 


(23) 


The results of this section will be used below to construct 
the center-of-mass dependence of the pairing amplitude 
in symmetric gauge. 


B. Andreev wavepacket scattering 


The dynamics of the problem we are interested in is 
generated by the Bogoliubov-de Gennes Hamiltonian op¬ 
erator 


T-LbcLG = 


(p- 


-E f 


At 


A 

(p+f a ) 2 


,(24) 


+ E f 


where the center-of-mass coordinate dependence and the 
relative coordinate of the pairing operator— can be ex¬ 
panded as 


<r'|A|r> = ]TA,*, 

3 

3 


r + r 


-j Xj(r~ r ') 


(25) 


d 2 k 


(27r) 


2 d i( k > 


3 ik-(r—r') 


(26) 


In what follows, the restriction in the sum over j will be 
made to the lowest term and, making use of Ea. E^l) . 

Ao^o(r) = A J2 e i7rril ” 2 e-5 c ACR e 2C^ e - 2 r€, ( 27 ) 

R 


For nodal d-wave superconductor 


do(k) 



(28) 



















For an s-wave superconductor the above quantity would 
be equal to unity. 

We are now in the position to define our scattering 
problem. The Nambu spinor \ip t ) evolves in time accord¬ 
ing to 





,h i 

= 'Hb<ig\iIh)- 

(29) 

We write 









HbcLG 

= Ho + V, 

(30) 

where 






Ho = ( 

huj c 

(at 

a + i) - 
0 

Ep 0 

— huj c [b^b + 

i) + B > 


0 

At 

<< <=> 

)■ 


(32) 

and separate 

the 

i time evolution due to Hq 

as 




Wt) = e 

~ iHot \m)- 

(33) 


Standard time dependent perturbation theory^ gives 


Wit)) 


i rt . . 

IV'(O)} + tt / dt'e* Hot Ve~^ Hot I Mt')) (34) 

^ Jo 

1 pt 

|^(0)) + — / dt'e^ Hot Ve~^ Hot IV’(O)) + ■ • ■ 

Jo 

(35) 


We are interested in finding (r| ij) t } given the initial state 
being the magnetic coherent state, which, without loss of 

I 


generality we choose to be purely hole-like 



This state is not an eigenstate of Hq, but its time evolu¬ 
tion due to Hq is known exactly. At t = 0, it corresponds 
to a Gaussian wavepacket peaked at (xo + iyo)/2£ = 
((3 q + iao)/y/2. The time evolution due to 7?o-only 
makes the complex variable «o time independent, and 
0o [t) = e lulct 0o- Therefore, ao determines the position of 
the center of the circle, and fio it) the radius of, and the 
angle along, the circle describing the classical motion of 
the wavepacket. The shape of the wavepacket does not 
change in time. 

We expect that when the initial wavepacket is prepared 
at an energy far away from the Fermi level, the effect 
of the pairing term is small, and that such wavepacket 
remains hole-like and that it continues to move along the 
circular trajectory. Therefore, if \huj c 0o0o — Ep\ ^ A, 
there should be no appreciable Andreev scattering. The 
goal is to determine the condition on ao and fio which 
would mark the transition from the regime where the 
wavepacket is unaffected by the condensate to the regime 
where the scattering is significant. Once such condition 
is identified, it will be utilized to define the energy scale 
e*, which is in turn used to achieve the scaling collapse 
of the non-perturbative numerical calculation for lattice 
d-wave superconductor in the mixed state. Although we 
are interested in the limit hoj c « A< Ep , this limit will 
be taken only at the end of the perturbative calculation. 


We imagine evolving the wavepacket from time 0 to time t, which is of order H/ A. Using the resolution of identity 
in terms of the magnetic coherent states (ED we find 


<«#t> 


0 i(ju c -£F/S)t 


0 

(r\a o ,e lulct 0o) 


— [ dt'e i< ' u>c ~ 2EF ^ h ^i t '~ t ) f 
ifr Jo J 


da\dai 

2iri 


dfitdfi i / (r|aie /5i|A|« 0 , e l “ ct '/3 0 ) 

2tti \ 0 


(37) 


The integral over a\ and /?i can be performed exactly, and so can the integral over k, the Fourier wavevector used to 
define ,\o( r ) in Eq. (E6l) . 


In order to perform the time integral in Eq. (1371) , we now take the limit 


flbj c „ Ep 


and assume that t < H/A. 


(38) 


After a somewhat lengthy but straightforward calculation (see Appendix for details), one finds that the dominant 









term for scattered part takes the form 


1 

ih 


dt ' e i^-2ErlW-t) I da^i I 


-{r\a 0 /3 0 (t)> ( X! e i7rnin2 e"5l^- 2 «l 2 e^C-CRr 
V R / 

j_ (Mt) m\ 


A 


2 ih \m A)(t) / -2t^ + 2iu> c PS(t)0o(t) 

where fio(t) = /3 0 e iu>ct . 


1 _ e ( 2i_ Tr -2*u; c 05(t)0o(t))A 


(39) 


Thus, in the stated limit, the wavepacket is appreciably 
Andreev scattered only if 


A 


/3 0 *M 


M 

fob) 


2 |J5 f - hu c P£(t)l3o(t)\ 


> 1. 


(40) 


Recall that the typical value of the kinetic momentum op¬ 
erator for the hole, p x +1 A x +i (p y + | Ay) , is \/2ihf3* /£. 
Therefore, we can interpret the term in the numera¬ 
tor in the above expression as the d-wave form factor 
amplitude. Such term is of course peaked in the anti- 
nodal direction. The term in the denominator repre¬ 
sents the difference between the typical energy of the hole 
wavepacket, i.e. the peak value of (p + -A) /2to, and 
the Fermi energy. The Andreev scattering is therefore 
maximized when the wavepacket is near the Fermi sur¬ 
face in the anti-nodal direction. If it is far from the Fermi 
surface, or is near the node, the wavepacket continues 
moving along the constant energy contours at A = 0, es¬ 
sentially as if the system was a normal metal. The above 


condition therefore marks the transition from the en¬ 
ergy regime where the wavepacket continues to move ac¬ 
cording to the semiclassical dynamics along the contours 
of constant energy, essentially undisturbed by the su¬ 
perconducting condensate, and the lower energy regime 
where the Andreev scattering occurs on time scales much 
shorter than 1 /u> c i.e. the time scale the wavepacket 
would need to complete the orbit. It is in the lower en¬ 
ergy regime that we find the suppression of the thermal 
Hall conductivity, as discussed in the next section. 


III. LATTICE D-WAVE NUMERICAL 
CALCULATION AND SCALING 

A. Tight-binding model 

In this section we resort to the numerical calculation of 
the intrinsic contribution to n xy along the lines discussed 
in Refs, d) and (0). 


We work on a two dimensional square lattice of spacing a - that we set to unity - and perpendicular magnetic field 
H. The tight-binding Hamiltonian describing the excitations is 


* = £ 
r 


( £ t r ,r+Scl t<T C r+ S,a + A r , r +,5 ^C^cJ + 5 ) j_ C r, 4 . C r+< 5 ,t) 

\< 5 =x,y 


H.C. 


- 



(41) 


Here c rjCT is the electron annihilation operator on the tight-binding lattice site r, not to be confused with the vortex 
lattice. 


The sum over the spin projection a =j“ or j, in the first 
and the last term of Eq. m is implicit; H.c. stands for 
Hermitian conjugation. The (nearest neightbor) hopping 
occurs in the presence of the uniform magnetic field, en¬ 
coded in the Peierls phase factor 


tr,r+s = (42) 


The magnetic flux $ through the elementary tight-biding 
plaquette appears through the link integral of the vector 


potential 


A r . r+ x = -7ry— , (43) 

00 

$ 

A r ,r+y = 7TX—— . (44) 

00 

The electronic flux quantum is 0o = hc/e. 

The ansatz for the tight-binding lattice pairing term is 


A rr j ^ 


A x 


A s e lB ^eiK +SdlV9 

-Ay = A, 


(45) 

(46) 


and the line integral is over the nearest neighbor link. 
Vortex positions, Rj, are inside the centers of some of 















the elementary plaquettes. They enter the pairing term 
through Q(r) which is chosen to be the solution of the 
continuum London’s equations 

V x V0(r) = 2 ttz ^2 <5(r - Rj) (47) 

j 

V • V0(r) = 0. (48) 

Vortices are positioned in the square lattice arrangement, 
with the vortex lattice at 45° relative to the underlying 
tight-binding lattice. As shown in Ref{H| the results dis¬ 
cussed below are largely independent of this choice. Each 
LxL magnetic unit cell is threaded by magnetic flux hc/e 
and contains a pair of vortices. Although the notation in 
this section uses the upper case letter L to denote the pe¬ 
riod of the magnetic unit cell, because the square vortex 
lattice is considered, in the tight-binding lattice units, it 
is equivalent to l B introduced earlier. 

The closed form solution of the London’s equations for 
the pairing field with such arrangement of vortices^, en¬ 
suring that the superfluid velocity, which is proportional 
to |V0(r) — -A(r), vanishes on average, is 

2 

0 0 ) = ( ar g [ CT ( 2 - ^;w,u/)] + 7^-2 ( z Zj - z * Zj )) . 

3 =1 

(49) 


In the above, just as in Ea. (l45l) . the line integral is again 

along the nearest neighbor link. The periodic factor 
(2) 

z r+8 r = 1 on each nearest neighbour link except the ones 

(2) 

intersecting the branch cut where 2 p A p = — 1. The iden¬ 
tity between the site factors on the left hand side of the 
Eq. (l5ll) and the link factors on the right hand side of (ICTl 
follows from considering products over the links forming 
closed clockwise loops around elementary tight-binding 
plaquettes. The left hand side must give +1 around 
each such elementary loop, regardless of whether such 
a loop contains a vortex, because it consists of a product 
of complex numbers with unit magnitude on each site. 
On the other hand, such a closed loop product formed 
from A dI ' ve mus ^ gj ve if fh e i o0 p contains a 
vortex and +1 if it does not, because j> dl ■ V0 = ±27t 
in the first case and j> d 1 ■ V0 = 0 in the second. For 
L = 28 considered here, inside the first Brilloin zone, 
the factor ei L di-ve can conve niently replaced by 

(1 + e iO(.r+d) e -iS(r)^ /| f + e iO(r+S) e ~i0(r) | _ Xhis wa ^ Qnly 

site variables enter the numerical calculation, and the 
link integrals need not be performed. 

The Heisenberg equations of motion 

ih-^1pr,a( k) = [V’r.o-Ck ),'H) 

= H BdG ( k)^ r , ff (k), (52) 


Here, 2 = x + iy (in tight-binding lattice units), Zj's 
denote the vortex positions inside the magnetic unit cell, 
and a(z; w, u/) is the Weierstrass a functions with periods 
oj = L and ui' = iL. 

The singular gauge transformation^’— ^- 4 turns the 
hopping and the pairing terms in Eg. (1411) periodic with 
the periodicity LxL, enabling the use of Bloch theorem. 
Performing the operator change of variables 


define the tight-binding lattice Bogoliubov-de^ Gennes 
single particle Bloch Hamiltonian operator, H B dG{ k), 
whose discrete eigenvalues, _E n (k), and eigenstates |nk), 
are labeled by the magnetic sub-band index n. For each 
k, there are 2L 2 such eigenstates. 

B. Thermal Hall conductivity 


c r ,f 

J 

C r,4 


— V(^ e ' 0 (r) ^ r ,t(k) ^ 

e-^W^(k) ) 


(50) 


where N uc is the number of magnetic unit cells in the 
entire lattice, ipr,o( k) is periodic in r with the periodicity 
of the magnetic unit cell, and k is within the magnetic 
Brilluoin zone ~x~ k x ,y < j;- 

The factors e i 9 ^' 1 must be handled with care due to 
the sign ambiguity associated with taking the square-root 
of a complex number. To start with, we connect vortices 
pairwise within each magnetic unit cell with branch-cuts, 
which are themselves periodic with the periodicity of the 
magnetic unit cell, and which intersect the elemenary 
tight-binding links. We chose the sign of the square-root 
such that the following identity holds 

e ie(H-«) e -* 0 (r) = ^ r ei i; +s di-ve_ (51) 


As mentioned at the end of the previous section, we 
denote by |nk) the eigenfunction of H B( ic{ k) with energy 
E n ( k) 

-ffsdG(k)|?rk) = _E„(k)|rak). (53) 

Then, the thermal Hall conductivity at temperature T 
has been shown to be given by^d£ 



where the Fermi occupation factor is 


e {/(fc B r) + 1 ’ 

and 


(55) 
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d 2 k 
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dii BdG (k) 
dk x 


?rk ) ( nk 


dH BdG (k) 

dk v 


mk \ — (x y) 


(k) 


(S m (k) - £„(k)) J 


(56) 


In the above, the double sum over the magnetic sub-band quantum labels m and n is to be performed subject to the 
stated restriction that for the given k, E m { k) < £ and E n ( k) > £. 


It is well known that the above formula can be written as 
the sum over occupied bands’ k-space integral over the 
Berry curvatur o 16 ’ 17 : 



where 


A„(k) = (nk|V k |nk). (58) 


For each fully occupied band, the integral extends over 
the entire magnetic Brillouin zone, and the occupied 
band contribution to C is an integer— d 6 ' 1 ' , the first Chern 
number. 

Therefore, determining the energy dependence of the 
k-space integral over the Berry curvature leads to finding 
the temperature dependence of the intrinsic thermal Hall 
conductivity. 

For the case of lattice d-wave superconductor con¬ 
sidered here, m, at H = 0 the pairing amplitude is 
Ak = 2A(cos k x — cos k y ) and the normal state disper¬ 
sion is ek = —2t(cosfca; + cos k y ). As discussed in the 
introduction, the condition 



£k ~ fJ- 


(59) 


results in two solutions (shown by red lines in FigJT]). Be¬ 
cause of the particle-hole asymmetry in the tight-binding 
dispersion, the value of |ek — y\ in the anti-nodal direc¬ 
tion along the two contours given by Ea. (l59l) is not the 
same. The lower of the two values of |ek — is 


e* 


A 4 ~- ImAI 

d* + |A/t| 


(60) 


The result of the numerical calculation for the model in 
Eci. m . with the temperature rescaled by e* with 0* = 
1 /v27t and the n xy with the value for A set to zero, is 
shown in Fig[2] Because the scaling with magnetic field 
has already been established, as was the independence on 
the vortex lattice geometry^, the above was computed for 
a single value of the magnetic length L = 28 and square 
vortex lattice. The method used here for an efficient 
computation of n xy has been detailed in RefjH. 


scale found in numerical calculations of the intrinsic ther¬ 
mal Hall conductivity in the mixed state of the nodal d- 
wave superconductor. Such picture is based the calcula¬ 
tion of the scattering of a magnetic coherent state within 
time-dependent perturbation theory and identifying an 
energy scale at which such scattering starts interfering 
with the ability of a wavepacket to complete its semi- 
classical orbit. Additionally, the results of the numerical 
calculation of K xy performed on a tight-binding lattice for 
the d-wave superconductor in the mixed state are shown 
to collapse well onto a single scaling curve (Fig©, pro¬ 
vided that the energy scale identified using the mentioned 
physical picture is used as the unit of temperature. These 
results show negligible dependence on the vortex core size 
as well as on the vortex lattice geometry. Such feature is 
also manifest within the wavepacket calculation. Similar 
calculation was performed in the case of a lattice s-wave 
superconductor, with on-site pairing term, which, unlike 
its d-wave counterpart, does not have k-dependence. In 
the s-wave case, the ^ dependence of the onset tempera¬ 
ture - which in the d-wave case amounted to 4 — |/r/t| - 
was absent. 

These findings may help establish measurements of n xy 
in very clean samples as a way to study the momentum 
structure of the pairing function in magnetic field via the 
bulk Hall transport method. 
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IV. SUMMARY 


The goal of this paper is to provide a physical picture 
which explains the existence of the onset temperature 















Appendix A: Details of the Andreev scattering of the magnetic coherent states 


A somewhat lengthy, but otherwise straightforward calculation, leads to 
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where 


( \ — IgjWo(t'-t) 
p = 2 ' 2 


1 + I e iaJ =(t'-t) i ’ 

J = i + igiMt'-t) {f*a> - 2C) 


J = 


ie 


iuj c (t'—t) 


1 + ie iw At'- 


— (y/le^po - 2r + Cr) • 


( 2 ) 

(3) 

(4) 


even when the entire (lengthy) expression is considered 
- is performed before the A integral. 

The above formula holds generally for any value of the 
ratio of the cyclotron frequency and the pairing ampli¬ 
tude. In order to perform the time integral in Ea. (|37l) of 
the main text, the limit of interest is taken 


hco c Ep 


(5) 


To obtain the above, first the overlap 
(«i, /3i|A|ao, e lLJct Po) is calculated in terms of the 
momentum integral; evaluation of the momentum 
integral is postponed until the the integrals over a\ and 
pi are performed. To calculate the integral over the 
momentum appearing in the d-wave form factor d( k), 
the denominator of (fc 2 — fc 2 )/k 2 is rewritten using the 

identity k~ 2 = dAe -Ak and the momentum integral, 
which is a product of a Gaussian and a polynomial 


It is also assumed that that the time duration does not 
exceed the time scale set by the pairing amplitude, i.e. 
that t < 1 /A. The terms containing complicated t depen¬ 
dence in the exponential can now be expanded to linear 
order in t. The resulting f-integrals are elementary. We 
postpone performing them for the sake of clarity, and in¬ 
stead rearrange the terms in order to reveal their physical 
content. Judiciously completing the squares, we find that 
the resulting expression can be brought into the form 
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( 12 ) 


where Po(t) = f3pe lu,ct and 


Similarly, 


Jo *(V2ia 0 -2 0 2 ’ 

J_ = _9_1_ (g) 

fo 4 ( v / 2/3oe i “= t -2e + CR) 2 ’ 

J -£ = ~l (V2ia 0 - 2?) (V2f3oe^ - 2£* + &) . (9) 


Eq.© has the form of a sum of two terms, representing 
the superposition of wavepackets. 

To analyze the first term in the Eq.©, note that the 
sum over R in the parenthesis corresponds to the super¬ 
position of Gaussians in £, whose centers are determined 
by the value of £r- The Gaussians are modulated by a 
pure phase factor. Therefore, if the value of £ is held 
fixed, then there is a value of R for which Cr comes close 
to maximizing the magnitude of the Gaussian. The next 
term multiplying the sum over R in the parenthesis is 
also a Gaussian in £ multiplied by a pure phase. It is 
peaked at 


Cpeafc — ^ 2 °° (I®) 

In the subsequent time integral, the values of £ and £* in 
the exponential are multiplied by a power of w c . Because 
the time interval is restricted to t < h/ A, the values of £ 
and £r may be replaced by their peak values inside the 
time integral 


>•peak 


- ('ipeak + iV 2a 0 ) = V2 ( ia 0 - /3q (t)) 


(ID 


1 1 _ 

Jo + Jo 2 ~ 8 + Po 2 (t) 

-t)[ia o 0o{t)-iy/2oi o C 

X e -^Ck(V2ao+2iC)uic(t'-t) ~ e iw c (t'-t)2/3*(t)/3 0 (t) ( 33 ) 

When the time integral is performed, the denominator 
containing Ep — Sw c /3 q /3o appears. In the state limit, 
this forces the entire expression to vanish, unless the 
value of ftp for the wavepacket of interest is such that 
hucPpfio ~ Ep. However, because of the term j q 2 + jg~ 2 , 
such scattered wavepacket is effectively suppressed by one 
power of hiOc/Ep. 

The second term Eq.© is peaked at 

Cpeafc = T ^/2^o(^) (I 4 ) 

which makes the sum over R dominated by the value 

Cr “* « 2£peafe = V 2 {iao + Po (t )), (15) 

allowing for the replacements 

11 9/1 1 \ 

ww + sw’ (6) 

— «J«(e)A(«), (17) 

P 0 U 

and 

e ~2iuj c (t'-t)(2C£+ia 0 @o(t)- V2(iaoC +iPo(t))) 
x e (k('/2°‘o+2iS)uc(t'-t) ~ e iu c (t'-t)2p£(t)p 0 (t) ^g) 

Because the factor jojo/po contains an additional factor 
of |/3o| 2 in the numerator, the suppression appearing in 
the first term discussed above is absent. Therefore, the 
dominant term is 
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( 19 ) 


This means that in the limit Say C A< Ep, over a time 
interval 0(H/A), the wavepacket is appreciably Andreev 


scattered only if 


Ep 
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■mmt) 


< 1 


( 20 ) 





















and if the d-wave form factor amplitude 


( 21 ) 


the constant (A = 0) energy contours, essentially as if the 
system was a normal metal. 


hit) flgm 
m) & (t)j 


is maximized i.e. the wavepacket is located in the anti¬ 
node. Otherwise, the wavepacket continues moving along 
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